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Abstract 

This  paper  describes  several  methods  of  extending  a  bivariate  interpolant  defined  on 
triangles  smoothly  beyond  the  triangulation.  Existing  methods  are  reviewed,  new  methods 
are  introduced,  and  the  methods  are  compared. 
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SIGNIFICANCE  AND  EXPLANATION 


The  general  problem  area  addressed  in  this  report  is  that  of  the  interpolation  of  scat¬ 
tered  bivariate  data  {xi,yi,Zi)  as  arise  for  example  in  the  Computer  Aided  Geometric 
Design  of  geometric  objects  like  the  shape  of  a  ship,  or  the  fuselage  and  wings  of  an  aircraft. 
Many  existing  techniques  proceed  by  triangulating  the  domain  and  then  defining  the  in- 
terpoiant  piecewise  on  each  triangle.  This  approaurh  has  the  drawback  that  the  interpolant 
cannot  be  readily  evaluated  outside  of  the  triangulated  domain.  This  report  deals  with 
the  problem  of  extending  such  a  triangular  interpolant.  Existing  methods  are  reviewed, 
new  methods  are  introduced,  and  the  methods  are  compared.  It  appears  that  at  present 
no  completely  satisfactory  techniques  exist  and  that  the  subject  requires  more  research. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive  summary 
lies  with  MRC,  and  not  with  the  author  of  this  report. 
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1.  Introduction 

The  interpolation  of  multivariate  scattered  data  is  a  rich  and  difficult  problem  area  of 
increasing  importance,  e.g.  in  the  Computer  Aided  Geometric  Design  of  geometric  objects 
like  the  shape  of  a  vehicle,  or  in  the  modeling  and  representation  of  measured  data.  For 
an  introduction  to  this  area  see  e.g.  the  survey  by  Barnhill,  1984,  the  proceedings  edited 
by  Barnhill  and  Nielson,  1984,  and  the  references  quoted  therein. 

In  this  paper,  we  restrict  ourselves  to  interpolation  of  scattered  data  in  two  indepen¬ 
dent  variables.  Thus  we  assume  that  we  are  given  points  {vertices)  V^,  =  (ij,y<)  €  31*,  t  = 
1 ,  •  •  • ,  iV  and  data  at  those  points.  The  data  are  either  part  of  the  problem  under  con¬ 
sideration,  or  else  they  have  to  be  made  up  suitably.  Most  frequently,  the  user  supplies 
only  function  {positional)  values.  Most  interpolation  schemes,  however,  also  require  values 
of  certain  derivatives  for  their  construction,  which  have  to  be  generated  from  the  given 
information. 

A  widely  used  approach  to  the  interpolation  problem  ( Triangular  Interpolation)  con¬ 
sists  of  first  triangulating  the  convex  hull  of  the  given  points,  D,  and  then  defining  the 
interpolant  piecewise  on  each  triangle,  taking  care  to  maintain  desired  global  properties 
such  ais  differentiability  of  a  required  order.  This  technique  has  several  advantages.  The 
interpolants  are  usually  of  a  simple  structure  (e.g.  polynomial),  and  they  can  be  defined 
and  evaluated  locally  (using  data  only  on  the  triangle  containing  the  point  of  evaluation). 
Much  of  the  analysis  and  construction  of  the  scheme  can  therefore  be  confined  to  a  sin¬ 
gle  triangle.  Many  triangular  interpolation  schemes  have  been  constructed,  see  the  above 
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references  for  details.  For  information  on  the  nontrivial  problem  of  constructing  a  trian¬ 
gulation  see  the  papers  by  Lawson,  1977,  Rarnhill  and  Little,  1984,  and  Cline  and  Renka, 
1984. 

A  drawback  of  triangular  interpolation,  however,  is  that  the  interpolant  cannot  be 
readily  evaluated  in  the  exterior,  E.  of  D.  The  notable  exception  to  this  general  statement 
is  provided  by  Akima’s  technique  (Akima.  1978),  which  is  very  efficent  and  widely  available, 
but  receives  only  a  “C”  grade  in  visual  quality  (Franke.  1978).  Akima  proceeds  by  dividing 
the  exterior  of  D  into  semi-infinite  rectangular  and  triangular  elements  and  then  defining  a 
special  structure  of  the  interpolant  piecewise  on  those  elements.  The  approach  is  described 
in  section  2.2  as  Discrete  Extension. 

In  the  same  survey,  Franke  extended  Nielson’s  minimum  norm  network  interpolant 
(Nielson,  1983)  beyond  £),  but  in  the  process  had  to  reduce  the  smoothness  of  the  inter¬ 
polant  (from  first  order  differentiability  to  mere  continuity).  Franke  uses  the  same  semi¬ 
infinite  elements  as  Akima,  and  extrapolates  by  Taylor  polynomials,  either  by  a  univariate 
one  perpendicularly  outward  from  edges  of  £>,  or  by  a  bivariate  one  about  a  boundary 
vertex  of  D.  This  technique  is  described  below  as  Transffnite  Extension.  Franke  did 
not  intend  his  extension  to  be  an  answer  to  the  triangular  extrapolation  problem,  rather 
he  had  to  modify  an  existing  technique  to  facilitate  comparison  with  others.  This  required 
only  very  little  extrapolation. 

In  this  paper,  we  describe  two  more  techniques  of  triangular  extrapolation.  One 
of  them.  Blending  by  Visible  Edges,  can  be  implemented  without  having  to  utilize 
the  particular  structure  of  the  interpolant.  All  that  is  required  is  the  ability  to  extend 
(smoothly)  the  representation  of  the  interpolant  on  any  particular  triangle  to  all  of 
This  is  possible  for  all  bivariate  triangular  interpoJants  known  to  this  author. 

The  second  new  technique.  Triangular  Extension,  consists  of  adding  vertices  to 
the  domain  structure,  redefining  the  triangulation,  generating  all  needed  data  at  the  new 
vertices,  and  then  redefining  the  interpolant  on  the  new  triangulation. 

Before  discussing  the  techniques  in  detail  we  mention  some  possible  issues  in  extending 
a  triangular  interpolant.  The  following  questions  may  be  relevant.  (4>  is  the  interpolant.) 

-1-  Is  the  degree  of  differentiability  of  $  on  E  at  least  the  same  as  that  on  D1 

-2-  Is  the  extended  domain  of  $  finite  or  does  it  cover  all  of  5R*? 

-3-  Is  0  of  the  same  structure  (e.g.  piecewise  polynomial  of  a  certain  degree)  on  E  as 
on  Z?? 

-4-  Is  the  class  of  functions  that  are  reproduced  exactly  by  the  scheme  maintained  in 
extending  from  D  to  E? 

-5-  Is  the  extended  version  of  identical  to  the  original  one  on  £)? 

-6-  Is  the  extension  technique  general  or  does  it  depend  on  the  special  structure  of 

on  Dl  For  example,  a  technique  is  not  general  if  it  requires  derivatives  of  since  these 
may  be  hard  to  supply  if  $  is  complicated  or  if  it  is  defined  by  a  black  box  routine. 

-7-  What  is  the  quality  of  the  technique  (e.g.  with  respect  to  accuracy  or  visual 
pleasantness)? 

The  following  table  lists  the  answers  to  the  first  six  of  the  above  questions  for  what  in 
this  author’s  opinion  constitutes  an  “ideal  technique”  and  the  three  techniques  described  in 
detail  below.  (Not  even  remotely  precise  criteria  exist  for  answering  the  seventh  question.) 


Question 

Ideal  Technique 

Transf.  Ext. 
(Franke) 

Discr.  Ext. 
(Akima) 

Blending  by 
vis.  Edges 

Triang.  Ext 

-1- 

yes 

no 

yes 

yes 

yes 

-2- 

all  of  5R2 

all  of  SR2 

all  of  5R2 

all  of 

finite 

-3- 

yes 

no 

no 

no 

yes 

-4- 

yes 

no 

no 

yes 

yes 

-5- 

yes 

yes 

yes 

yes 

depends 

-6- 

general 

special 

special 

general 

special 

Table  1:  Properties  of  Methods 


It  is  apparent  that  none  of  our  techniques  is  “ideal’'.  The  choice  of  an  extrapolation 
technique  depends  on  the  problem  under  investigation  and  the  underlying  interpolant  on 
D.  Indeed,  if  evaluation  outside  of  D  is  necessary,  then  an  alternative  method  that  does  not 
require  a  triangulation  (e.g.  one  of  those  described  in  Franke,  1982)  may  be  preferable  in 
spite  of  the  advantages  and  attractive  features  of  triangular  interpolation.  Note,  however, 
that  the  mere  ability  to  evaluate  an  interpolant  outside  of  D  does  not  by  itself  make  that 
method  a  worthy  extrapolation  technique. 

8.  Extrapolation  Techniques 
8.0  Geometry  and  Notation 

Before  entering  the  discussion  of  the  individual  extrapolation  techniques  we  need  to 
introduce  some  additional  notation.  We  assume  that  D  is  the  convex  hull  of  the  given 
vertices  K,,  i  =  and  that  it  has  been  triangulated.  Of  particular  interest  are 

boundary  vertices  of  D  and  triangles  having  at  least  one  edge  on  the  boundary  of  D. 
Without  restriction  of  generality  we  assume  that  the  vertices  have  been  labeled  such  that 
V,-,  i  =  are  the  boundary  vertices  in  counterclockwise  order.  The  boundary 

segment  with  endpoints  Vi  and  V,+i  (where  :=  V'])  is  denoted  by  e,.  The  unique 

triangle  of  which  «,  is  an  edge  is  Aj.  (If  a  triangle  has  two  boundary  edges,  e,  and  e,+i, 
then  Ai  =  A,+  j.) 

We  denote  the  extended  interpolant  by  On  a  triangle  A,,  it  is  represented  by 
♦i.  We  assume  that  ♦i  can  be  evaluated  anywhere  in  SR^.  This  stipulation  is  necessarily 
somewhat  vague  because  of  the  abundance  and  variety  of  available  interpolation  schemes. 
For  any  given  method,  it  will  usually  be  obvious  how  to  proceed.  For  example,  if  a  scheme 
is  polynomial  on  each  triangle,  then  will  be  the  same  polynomial  in  all  of  3?^. 

The  geometric  concepts  introduced  in  this  paragraph  are  illustrated  in  figure  1.  The 
exterior  of  D  is  E.  We  divide  E  into  semi-infinite  regions  by  drawing  two  lines  from  each 
boundary  vertex  that  are  perpendicular  to  the  adjoining  two  boundary  edges.  The  semi¬ 
infinite  rectangular  element  corresponding  to  c,  is  7?,.  and  the  wedge  corresponding  to 
is  W,.  A  wedge  W,-  is  empty  if  the  two  boundary  segments  joining  at  V,  are  parallel.  The 
triangles  in  D  and  the  set  of  Ri  and  W,  tesselate  the  plane. 
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For  a  given  point  P  e  we  denote  by  B  the  perpendicular  projection  of  P  onto  Cj. 
We  describe  the  location  of  P  in  terms  of  new  coordinates  s  and  t  as 


P  =  B- 


HP-B) 

\\P-  B\ 


where  B  is  the  convex  combination 

B  =  sV.+i  -  (1  -  s)V,  =  V,  +  s(V.^,  -  V.) 
and  s  and  t  are  given  explicitly  by 

^  ^  =  B\\  (1) 

In  any  implementation  of  a  scheme  based  on  the  above  tesselation  one  must  be  able 
to  compute  in  which  element  a  given  point  resides.  This  can  be  accomplished  as  follows: 

Location  Algorithm 

Note:  Ail  index  arithmetic  is  modulo{M). 

1.  Given  a  point  P  €  Determine  a  specific  triangle  in  £>,  call  it  the  current  triangle. 

2.  Compute  the  barycentric  coordinates  bi,  f>2>  and  63  of  P  with  respect  to  the  current 
triangle  (barycentric  coordinates  are  discussed  in  several  of  the  articles  in  Barnhill 
and  Nielson,  1984). 

3.  IF  all  three  barycentric  coordinates  are  non-negative,  THEN  P  resides  in  the  cur¬ 
rent  triangle,  so  STOP,  ELSE  consider  that  edge  e  of  the  current  triangle  that 
corresponds  to  the  most  negative  barycentric  coordinate.  //’  e  =  is  a  boundary 
edge,  THEN  P  lies  outside  of  D,  GO  TO  step  4.  IF  e  is  an  interior  edge,  THEN 
replace  the  current  triangle  by  the  triangle  across  e  and  GO  TO  step  2. 

4.  Let  Sold  =  1/2.  GO  TO  step  5. 

5.  Using  (1)  compute  s.  IF  s  e  [0,  IJ  THEN  P  €  Rx  hence  STOP.  IF  s  <  0  and 
Sold  >  1  then  P  G  W,_i  hence  STOP.  IF  s  >  1  and  Soid  <  0  THEN  P  € 
hence  STOP.  IF  none  of  the  above,  THEN  GO  TO  step  6. 

6.  Let  Sold  =  'S.  IF  s  <  0  THEN  i  :=  t  -  I  ELSE  i  :=  i  -i-  1  .  GO  TO  step  5. 

2.1  Transfinite  Extension 

The  general  idea  of  Franke’s  extension  is  ais  follows.  In  each  wedge  Wj  define  $  by 
bivariate  Taylor  interpolation  of  order  m,  say,  to  the  derivatives  of  $  (as  they  are  defined  on 
D)  at  the  vertex  V,.  In  each  rectangle  i2,  define  $  by  Taylor  extrapolation  to  derivatives  of 
$  along  lines  perpendicular  to  the  boundary  edge.  Note  that  it  must  be  possible  to  evaluate 
derivatives  of  $  at  any  point  on  the  boundary  (hence  the  name  transfinite  extension).  The 
degree  of  the  transfinite  extrapolation  must  equal  m,  for  otherwise  the  extended  interpolant 
would  be  discontinuous  along  rays  joining  a  wedge  and  a  rectangle.  Also,  the  derivative 
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data  must  be  well  defined,  i.e.  ^  must  be  at  least  m  times  differentiable  at  the  boundary 
vertices.  (We  assume  that  it  is  arbitrarily  often  differentiable  across  a  boundary  edge,  so 
that  all  needed  derivatives  can  be  supplied.) 

To  be  explicit,  we  list  the  appropriate  formulas.  Let  P  =  (x,y)  €  SR*  be  a  general 
point.  If  P  e  Wi  we  define 


♦(/>):=  y; 


li+u<m 


(Vi)!!  -  -  y,)‘ 


(2) 


If  P  €  Ri  we  express  the  location  of  P  in  terms  of  the  new  coordinates  t  and  s  as 
described  in  (1)  and  define 

=  (3) 


where  partial  differentiation  denotes  the  standard  normalized  directional  derivative. 

The  interpolant  extended  in  this  fashion  will  be  m  times  differentiable  everywhere  on 
the  boundary  of  D.  However,  the  smoothness  of  $  may  be  reduced  across  the  rays  sepa¬ 
rating  semi-infinite  elements.  Assume  that  P  lies  on  a  ray  emanating  from  the  boundary 
vertex  V,  and  that  a  ib-th  order  partial  derivative  is  to  be  evaluated  at  P.  Differenti¬ 
ating  in  (2)  and  (3)  yields  derivatives  of  ^  through  {m  +  k)  -  th  order  at  B  =  Vi.  Thus, 
in  order  for  to  be  continuous,  ^  must  be  at  least  (m  -h  k)  times  differentiable  at  Vi. 

Consider  for  example  Akima’s  interpolamt,  which  is  twice  differentiable  at  the  vertices. 
There  one  could  pick  m  =  1,  i.e.  interpolate  linearly  in  triangular  wedges,  and  linearly 
along  rays  emanating  perpendicularly  from  the  boundary  segements.  The  resulting  ex¬ 
tended  interpolant  would  then  be  once  differentiable  everywhere  in  SH*.  This  property, 
however,  is  crucially  dependent  upon  the  fact  that  Akima’s  interpolant,  while  only  once 
differentiable  on  edges  separating  triangles,  is  actually  twice  differentiable  at  the  vertices 
of  the  triangulation.  Akima  does  extend  his  interpolant  differentiably  to  all  of  91*,  but  he 
proceeds  in  a  different  manner  which  is  described  in  the  next  subsection. 


S.S  Discrete  Extension 

The  basic  idea  of  this  technique  is  to  interpolate  to  discrete  data  at  the  boundary 
vertices  and  to  construct  a  special  extension  on  each  of  the  two  types  of  semi-infinite 
elements.  The  vertex  data  are  read  off  the  interpolant  as  it  is  defined  on  D.  We  illustrate 
the  technique  by  considering  Akima’s  interpolant. 

In  a  triangular  wedge  W,  Akima  simply  uses  bivariate  quadratic  Taylor  interpolation  as 
defined  in  (2).  In  a  rectangular  strip  Ri,  however,  he  interpolates  by  the  unique  polynomial 
in  span{l,  s.  s*,  s*,  s*,  a®,  t,  st,  s^t,  a*<,  t*,  and  3a*t*  -  25*t*}  that  interpolates  to  the 
twelve  data  given  at  V,  and  V,+  i. 

Akima's  extended  interpolant  is  continous  across  e,  because  it  can  be  represented 
there  as  a  univariate  quintic  that  is  determined  uniquely  by  the  six  tangential  data  at 
V,  and  Similarly,  it  is  continuous  across  the  rays  separating  semi-infinite  regions 

because  there  it  can  be  represented  uniquely  as  a  univariate  quadratic  function  that  is 


i 


I 


».  J 


t 


-5- 


t 


I 


determined  by  the  vertex  data.  Differentiability  follows  because  a  perpendicular  cross- 
boundary  derivative  is  a  unique  univariate  cubic  across  boundary  edges,  and  a  univariate 
linear  function  across  rays. 

Note  that  the  construction  of  Akima’s  interpolant  depends  upon  the  special  structure 
of  $  on  D.  There  is  no  obvious  way  of  generalizing  this  approach. 


2.S  Blending  by  Visible  Edges 

Here  we  proceed  as  follows:  For  a  given  point  P  G  E  we  consider  all  visible  edges  Cj, 
i.e.  all  edges  for  which  the  determinant  of  the  matrix 


( 


1 

V. 


1 

P 


is  p(»itive.  The  weight  attached  to  a  triangle  A,  is  larger  the  larger  the  angle  a,  which  is 
formed  by  V^,  P  and  (see  figure  1).  The  weights  are  chosen  to  force  a  certain  degree 
of  global  smoothness  and  to  maintain  the  degree  of  precision  of  the  original  interpolant. 
We  define 

t 

where  the  summation  is  taken  over  all  visible  edges  and  the  blending  functions  'J'(at) 
satisfy 

£4-(a.(P))  =  1.  (S) 

i 

Obviously,  if  $i(P)  =  F(P)  for  some  given  function  F  and  all  i,  then  ^{P)  =  F{P) 
so  that  the  precision  of  the  interpolant  on  D  is  preserved. 

We  next  address  the  question  of  the  smoothness  of  $  as  defined  in  (4).  Let  us  assume 
that  the  individual  terms  in  (4)  are  m  times  differentiable.  Then,  applying  any  m-th  order 
differentiation  operator  to  $  yields  a  sum  of  products  of  derivatives  through  order  m  of 
’J',  and  ttj,  all  of  which  are  continuous  by  assumption.  However,  if  some  of  the  angles 
a,(P)  approach  zero,  then  the  corresponding  term  will  be  added  or  omited  from  equation 

(4) .  In  order  to  maintain  continuity  of  the  derivative,  all  corresponding  terms  must  be 
zero  when  this  happens.  Each  term  contains  some  derivative  of  ^  as  a  factor.  Thus  we 
require  that 

=  0  Vi' =  0,l,...,m.  (6) 

Note  that  by  picking  m  larger  than  the  degree  of  smoothness  of  $  on  £?  it  is  possible  to 
achieve  a  degree  of  smoothness  in  the  exterior  E  that  exceeds  that  in  D\ 

The  task  now  remains  of  constructing  blending  functions  that  satisfy  the  requirements 

(5)  and  (6).  The  simplest  approach  is  a  s  follows.  We  force  (6)  by  defining  functions 


X(a)  = 
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and  then  define  the  blending  functions  'ifi  by  the  standard  trick  of  dividing  by  the  sum  of 
the  appropriate  terms,  i.e. 


^i{P) 


x(c^i{P)) 

HjXiajiP)) 


(7) 


where  the  sum  in  the  denominator  is  taken  over  all  visible  edges. 


2.4  Triangular  Extension 

The  general  philosophy  is  as  follows;  Virtually  all  triangular  interpolants  require  for 
their  construction  certain  auxiliary  data  that  are  not  supplied  by  the  user.  Hence  there 
must  be  a  method  of  generating  such  information  from  user  supplied  data.  So  extend 
the  triangulated  domain  by  adding  points  in  the  domain  and  generate  the  data  needed 
at  the  additional  points  (including  the  kind  usually  supplied  by  the  user)  by  using  the 
data  generation  method  that  is  already  in  existence.  The  triangulation  of  the  extended 
domain  can  be  afcomplished  by  adding  triangles  to  the  existing  triangulation  of  D,  or  by 
triangulating  the  convex  hull  of  the  new  extended  point  set. 

Of  course,  not  every  method  for  the  generation  of  auxiliary  data  can  actually  be 
generalized  to  work  on  the  expanded  point  set.  Most  frequently,  the  user  supplies  C® 
data  and  the  interpolation  scheme  must  then  make  up  some  kind  of  derivative  data.  This 
process  amounts  to  numerical  differentiation,  and  is  clearly  not  suitable  or  generalizable  to 
the  generation  of  additional  positional  data.  However,  here  we  describe  a  specific  approach 
that  generalizes  easily  and  naturally. 

In  Alfcld,  1984,  as  well  as  in  earlier  references  (Schmidt,  1982;  Nielson,  1983)  the 
interpolation  problem  was  solved  by  picking  that  function  in  a  certain  space  of  functions 
that  minimized  a  particular  functional  subject  to  interpolation.  The  fundamental  idea  is 
to  consider  needed  data  (such  as  derivative  values  at  given  vertices  or  positional  values  at 
additional  vertices)  to  be  parameters  that  enter  the  minimization. 

The  details  of  the  algorithm  used  for  the  numerical  examples  in  section  3.  will  be 
reported  elsewhere.  However,  we  sketch  its  main  ingredients.  The  algorithm  is  applicable 
for  any  triangulated  domain  where,  at  any  vertex  V,  some,  none,  or  all  of  the  following 
data  may  be  given:  {/"(V),  F*(V),  Fy(V)}.  Since  it  is  possible  to  have  no  data  at  a 
given  vertex,  it  is  simple  to  extend  the  given  domain  and  define  the  interpolant  on  the 
retriangulated  extended  domain.  The  new  domain  is  of  the  same  structure  as  the  old  one, 
and  the  distinction  between  D  and  E  no  longer  applies. 

The  interpolant  is  cubic  on  each  triangle  (as  in  Schmidt,  1982),  and  is  globally  once 
differentiable.  The  number  of  data  supplied  must  not  exceed  the  total  number  of  available 
parameters.  Among  all  interpolants  q,  the  one  is  chosen  that  minimizes  the  thin  plate 
functional 

^(9)  =  <ilz  +  2qly  +  qlydxdy 


The  minimizing  solution  is  found  by  applying  the  finite  element  technique  of  setting 
up  element  stiffness  matrices  and  assembling  them  into  a  global  stiffness  system.  The 
interpolant  has  the  advantage  of  greatest  possible  simplicity  (i.e.  piecewise  cubic)  but 


has  the  drawback  that  the  differentiability  conditions  cannot  be  imposed  locally  (as  they 
could  in  Alfeld,  1984)  but  rather  have  to  be  incorporated  as  global  constraints.  It  is  well 
known  that  the  C*  conditions  may  be  redundant  (Schumaker,  1984,  Morgan  and  Scott, 
1977).  The  only  known  source  of  redundancy  are  singular  vertices,  i.e.  points  at  which 
exactly  four  triangles  meet  at  the  intersection  of  the  two  diagonals  of  the  quadrilateral 
defined  by  the  four  triangles.  If  a  singular  vertex  (or  nearly  singular  vertex)  occurs  then 
the  singularity  can  be  destroyed  by  introducing  an  additional  vertex  (without  any  data 
attached  to  it)  nearby  and  considering  the  data  needed  there  as  three  more  parameters 
that  enter  the  minimization  problem. 

The  resulting  linear  system  is  solved  directly,  rather  than  iteratively.  Powerful  meth* 
ods  for  this  task  exist  (George  and  Liu,  1981). 

S.  Numerical  Examples 

In  this  section,  we  illustrate  the  extrapolation  techniques  described  in  the  preceding 
section  by  interpolating  to  Franke’s  well-known  test  function 

/(i,y)  =  >/< 

4 

^?e-(9*+l)V49-(9v+l)/10 

4 

+  l--((9*-7p  +  (9v-3)=)/4 

_1  (9*-4)»-(9y-7)» 

5 

The  domain  data  and  the  triangulation  of  D  are  given  in  table  1.  At  each  vertex, 
the  function  value,  and  no  derivative  values,  are  given.  The  domain  D  is  the  unit  square 
covered  by  36  (essentially  random)  points  and  54  triangles.  The  extended  domain  was 
taken  to  be  a  square  with  twice  the  area  of  £>,  with  D  in  its  center. 

There  are  five  groups  of  plots  (which  were  generated  by  (PLOT79),  see  Beebe,  1979). 

1.  Primitive  Function 

2.  IMSL  Interpolant  This  interpolant  was  generated  by  the  IMSL  routine  IQH- 
SCV,  which  is  an  implementation  of  Akima’s  technique. 

3.  Blending  by  Visible  Edges  This  is  the  result  of  the  technique  described  in 
section  2.3. 

4.  36-1-4  point  extension  Here  the  domain  D  was  extended  by  adding  the  four 
vertices  of  the  extended  domain  E,  and  retriangulating.  Then  the  technique  described 
in  section  2.4  was  applied.  In  this  example  it  so  happened  that  the  triangulation  of 
the  original  domain  D  was  preserved,  although  this  cannot  be  expected  to  be  true  in 
general. 

5.  36  -f  12  point  extension  Similar  to  4.,  except  that  12  points  equally  spaced  on 
the  boundary  of  E  were  added.  As  in  4.,  the  triangulation  of  D  was  preserved. 

For  each  of  the  above  groups  there  are  two  plots,  namely  a  contour  plot  with  the  tri¬ 
angulation  superimposed,  and  a  hidden  line  view  of  the  surface  from  a  Northeast  direction. 
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The  plotting  parameters  (viewpoint,  vertical  scaling,  interval  between  contour  lines,  and 
range  covered  by  the  contour  lines]  are  identical  in  all  five  groups. 

It  is  apparent  that  none  of  the  extrapolated  interpolants  matches  the  simplicity  and 
beauty  of  the  primitive  function.  The  individual  surfaces  differ  substantially,  both  quan¬ 
titatively  and  qualitatively.  Blending  by  Visible  Edges  seems  to  give  rise  to  the  worst 
undulations,  introducing  several  spurious  saddle  points.  This  is  unfortunate  because  it 
offsets  the  other  advantages  of  this  method,  namely  generality,  versatility,  and  the  ability 
to  generate  any  degree  of  smoothness  in  E  (by  choosing  the  exponent  in  (7)  suitably  large). 
It  is  an  open  question  whether  better  results  can  be  obtained  by  a  different  choice  of  the 
blending  functions.  Akima’s  method  performs  better,  but  it  suffers  from  poor  quality  even 
in  D  due  to  deficiencies  in  the  derivative  generation.  Triangular  Extension  gives  rise  to 
surfaces  with  few  undulations  and  seems  to  be  the  best  of  all  methods.  The  quality  of  the 
surface  seems  to  improve  with  the  number  of  points  added  to  the  original  set. 


Primitive  Function 


Northeast  View 


IMSL  interpolant 


IMSL  interpolant 


Northeast  View 


36  +  4  point  extension 


Northeast  View 
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36  +  12  point  extension 


Conclusions 


Triangular  Extrapolation,  as  well  as  multivariate  extrapolation  b2ised  on  interpolants 
other  than  triangular  ones,  is  a  difficult  area  that  requires  more  research.  At  present, 
several  techniques  exist,  but  none  of  them  is  entirely  satisfactory.  Of  the  existing  ones, 
Triangular  Extension  seems  to  generate  the  visually  most  pleasing  surfaces. 
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